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ABSTRACT 


ASC  DEQ4  is  a  floating  point  subroutine,  written  in  FORTRAN  IV 
source  language,  which  integrates  numerically  a  set  of  N  simultaneous 
second-order  ordinary  differential  equations  in  which  first  derivatives  may 
or  may  not  appear  (i.e.  ,  y"  -  f(t,  y^,  yf  )  of  y'/  =  f(t,  y.),  i  =  1 ,  2,  ....  Nj. 

If  the  N  equations  can  be  separated  into  two  groups  (IB  and  N-IB)  such 
that  the  first  IB  equations  are  not  dependent  on  the  final  N-IB  equations 
(e.  g.  ,  variational  equations)  then  DEQ4  has  the  capability  of  integrating 
the  final  N-IB  equations  at  a  larger  step  size  than  the  first  IB  equations, 
thus  saving  2(R-1)(N-IB)  derivatives  per  integration  step.  This  subroutine 
obsoletes  subroutine  DEQ2  with  the  following  improvements:  better  accuracy 
controls,  new  starting  procedure,  improved  halving  and  doubling  procedure, 
reduction  in  computing  time,  and  reduction  in  core  storage  requirements 
(ION  less). 

The  subroutine  is  restricted  in  that  it  contains  20  digit  octal  constants 
(real  constants)  for  the  CDC  o000  series  machines. 
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SECTION  I 


INTRODUCTION 


A.  PURPOSE 

The  purpose  of  ASC  DEQ4  is  to  integrate  numerically  a  set  of  N 

simultaneous  second-order  ordinary  differential  equations  in  which  first 

derivatives  may  or  may  not  appear  [i.  e.  ,  yV  =  f(t,  y.,  y!)  or 

y'.1  =  f(t,  y.)  is  1,  2,  N].  If  the  N  equations  can  be  separated  into  two 

groups  (IB  and  N-IB)  such  that  the  first  IB  equations  are  not  dependent  on 

the  final  N-IB  equations,  then  DEQ4  has  the  capability  of  integrating  the 

K 

final  N-IB  equations  at  a  larger  step  size  (ratio  =  R  =  2  ,  K  £  1)  than  the 
first  IB  equations.  This  will  save  2(R-1)(N-IB)  derivatives  per  integration 
step  when  operating  in  the  Cowell  mode. 

B.  RESTRICTIONS 


The  restrictions  associated  with  the  use  of  ASC  DEQ4  are  as  follows: 

1.  No  internal  checks  are  made  for  overflow  or  underflow. 

2.  The  user  must  provide  an  auxiliary  subroutine  (DAUX)  which 
evaluates  the  second-order  derivatives  and  the  name  of  the  aux¬ 
iliary  subroutine  must  be  defined  as  an  argument  in  the  calling 
sequence  and  also  through  the  use  of  the  EXTERNAL  statement 
in  the  main  program  (i.  e.  ,  EXTERNAL  DAUX). 

3.  If  the  two-group,  multi-step  mode  is  used,  then  DAUX  should 
be  programmed  to  skip  the  evaluation  of  the  final  N-IB 
derivatives  when  TEST  is  negative.  (See  IB,  TEST,  and 
YMIN  in  Section  II- A  for  restrictions  on  the  use  of  the  two- 
group,  multi-step  mode. 

4.  Initial  conditions  must  be  stored  prior  to  entering  the  program. 

5.  ASC  DEQ4  is  a  single  precision  subroutine  written  in  Fortran 
IY  Source  Language. 

6.  The  subroutine  contains  20  digit  octal  constants  (real  constants) 
for  the  CDC  6000  series  machines. 


c. 


METHOD 


A  fourth-order  Runge-Kutta  method*  is  used  to  start  the  integration  and 
to  halve  the  step  size  during  integration.  A  Cowell  "second-sum"  method 
based  on  eighth  differences  is  used  to  continue  the  integration.  Doubling  is 
accomplished  in  the  Cowell  mode  through  the  accumulation  of  alternate  steps. 
Truncation  error  can  be  controlled  by  choosing  an  appropriate  step  size, 
and  by  using  the  variable  step  size  mode  of  operation.  DEQ4  has  an  automatic 
restart  procedure  which  works  as  follows: 

After  eight  Runge-Kutta  steps  (h  =  H/IR),  an  attempt  is 
made  to  integrate  in  the  Cowell  mode  (h  =  H/IR).  If  the 
error  criterion  (see  ER  in  Section  II-A)  cannot  be  sat¬ 
isfied  for  the  first  IB  equations,  the  initial  conditions 
(i.  e.  ,  tQ,  y£0>  and  yj^)  are  restored,  the  step  size  is 
reduced  to  h  =  h/2*IR),  and  the  integration  is  restarted. 

This  procedure  continues  until  the  error  criterion  is 
satisfied,  or  until  h  <  HMIN.  For  the  latter,  the  sub¬ 
routine  sets  TEST  =  13  and  exits  to  the  user. 

Note:  Either  12  <  0,  or  HMIN  =  HMAX  nullifies  the 
restart  procedure  as  well  as  the  halving  procedure. 


*J.  B.  Scarborough,  Numerical  Mathematical  Analysis,  Third  Edition, 
John  Hopkins  Press,  Baltimore  (1^55)  pp.  301-302. 
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SECTION  II 


PROGRAM  USAGE 


A.  CALLING  SEQUENCE1 

CALL  DEQ4(N,  II,  12,  LA,  IB,  IR,  ER,  HMIN,  HMAX,  YMIN, 

DAUX,  TEST,  IDH,  NTRY,  JHH,  JHD,  VMIN,  VMAX,  T, 

H,  Y,  YP,  Y2P,  Ti,  T2,  T3,  T4,  T5,  T6,  T7,  F2P,  FiP, 

DLT 1,  DLT2,  DLT3,  DLT4.  DLT5,  DLT6,  DLT7,  DLT8) 

The  nomenclature  is  as  follows: 

N  the  number  of  equations. 

11  an  option,  so  that  if 

II  >  0,  ist  derivatives  are  present  in  the  evaluation  of 
the  second  derivatives,  [yl1  =  f(t,  y^,  yp]> 

11  <  0,  ist  derivatives  are  missing  in  the  evaluation  of 

the  second  derivatives,  [yl1  =  f(t,  y^)]. 

12  an  option,  so  that  if 

12  £  0,  variable  step  size  mode  of  operation  is  used. 

12  <  0,  fixed  step  size  mode  of  operation  is  used. 

(Note:  For  12  <  0,  h  =  H/IR  for  all  steps.  ) 

IA  an  indicator  switch  to  the  user  during  exit  to  DAUX. 

LA  =  -1  for  first  '  1 ,  2,  3  if  Runge-Kutta)  pass  through  DAUX. 

IA  =  +i  for  final  (4th  if  Runge-Kutta)  pass  through  DAUX. 

This  applies  for  each  integration  step.  In  the  Cowell  mode, 
IA  =  -1  when  the  derivatives  of  the  predictor  are  being 
asked  for,  and  LA  =  +1  when  the  derivatives  of  the  corrector 
are  being  asked  for. 

IB  Only  the  first  IB  (£  N)  equations  are  tested  to  determine  whether 

it  is  necessary  to  halve  or  possible  to  double  the  step  size  or  to 
proceed  with  a  Cowell  integration  step.  (See  Section  I-A  and 
Section  II-C  for  an  additional  use  of  IB.  ) 

See  NTRY  and  Section  II-B 
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IR 


ER 


HMIN 


HMAX 


YMIN 


DAUX 


For  a  given  step-size  H,  the  initial  step  size  for  both  Runge- 
Kutta  and  Cowell  will  be  h  =  H/IR.  If  halving  is  required,  the 
current  step-size  h  is  reduced  to  h  =  h/(2*IR)  for  all  N  equations, 
and  the  integration  procedure  returns  to  the  Runge-Kutta  mode. 
(Usually  IR  =  8  or  16;  if  IR  =  0,  DEQ4  sets  IR  =  16.  ) 

The  user  should  set  ER  =  IE  -  S  (i.  e.  ,  10“^),  where  S  is  the 
approximate  number  of  significant  figures  desired.  If  ER  =  0, 
DEQ4  sets  ER  =  IE  -  11.  (User  should  test  for  best  ER.  ) 

In  general,  ER  £  IE  -  11  is  recommended. 

the  minimum  absolute  value  of  the  step  size  allowed  when 
halving  is  required.  During  the  starting  procedure  (see 
Section  I-C),  if  halving  is  required  and  h/(2*IR)  <  HMIN, 
the  integration  is  terminated  with  TEST  =  13.;  otherwise, 
the  integration  is  continued  in  the  Cowell  mode  with  the 
current  step-size  h  (i.  e.  ,  halving  is  not  permitted).  If 
HMIN  =  HMAX,  all  halving  is  suppressed  during  the  entire 
integration  procedure.  If  HMIN  =  0,  DEQ4  sets  HMIN  =  IE  -  5. 
Only  the  first  IB  equations  are  tested  for  halving. 

the  maximum  absolute  value  of  the  step  size  allowed  for  the 
first  J.B  equations.  The  final  N-IB  equations  may  be  integrated 
at  a  larger  step  size  when  operating  in  the  two-group,  multi- 
step  integration  mode  (see  Section  I-A).  If  HMAX  =  0,  DEQ4 
sets  HMAX  =  1.  A  value  of  HMAX  £  16.  is  recommended  for 
most  problems. 

the  minimum  absolute  value  allowed  for  y.  for  the  relative  error 
test  during  the  halving  and  doubling  procedure.  YMIN  prevents 
unnecessary  halving  for  y-  =:  0.  If  YMIN  =  0,  DECK  sets 
YMIN  =  1.  (User  should  supply  YMIN). 

Note:  For  the  more  difficult  problems,  it  may  be 
necessary  for  the  user  to  modify  the  subroutine  and 
make  YMIN  a  vector  of  dimension  IB  --  thus  allowing 
a  different  YMIN  for  the  first  IB  equations. 

the  location  of  the  entry  point  of  an  external  subroutine  (supplied 
by  the  user)  which  evaluates  and  stores  (see  Y2P)  the  second- 
order  derivatives  yt'.  DAUX  must  be  defined  by  the  user  through 
the  use  of  the  EXTERNAL  statement  in  the  main  program,  and 
COMMON  must  be  used  as  a  means  of  data  linkage.  CALL 
DAUX  is  used  by  DEQ4  to  enter  DAUX.  If  the  two-group, 
multi-step  mode  is  used,  then  DAUX  should  be  programmed  to 
skip  the  evaluation  of  the  final  N-IB  derivatives  when  TEST  is 
negative. 
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TEST 


has  the  following  multiple  uses: 

Initially  (i.  e.  ,  prior  to  NTRY  =  1),  the  user  must  set: 

TEST  =  +1.,  to  integrate  the  N  equations  at  the  same 
step  size. 

TEST  2  +2.,  to  use  the  two-group,  multi-step  mode.  The 
maximum  ratio  of  the  step  sizes  for  the  two  groups  will 
be  R-MAX  =  2 **( TES T  -1).  Recommend  TEST  £  5. 

After  each  integration  step  (NTRY  -  2)  DEQ4  will  set: 

TEST  =  +1.,  if  the  integration  was  a  Runge-Kutta  step 
(i.  e.  ,  during  starting  procedure,  or  halving  procedure 
if  JHH  =  1). 

TEST  =  +2.,  if  the  integration  has  been  restarted  during 
the  initial  starting  procedure.  Indicates  h  =  h/(2*IR)  and 
Runge-Kutta  step. 

TEST  =  +3.,  if  the  step  size  has  been  reduced  [h  =  h  /  ( 2 *IR , ] 
and  the  integration  has  been  returned  to  a  previous  step 
during  the  halving  procedure  (possible  only  if  JHH  =  3). 

TEST  =  -1.,  if  the  integration  was  a  Cowell  step. 

During  the  starting  procedure  (see  Section  I-C): 

TEST  =  +13.,  if  the  integration  has  been  terminated  during 
the  restarting  procedure.  Indicates  programming  error, 
or  ER,  H,  or  HMIN  too  restrictive. 

During  transfers  to  DAUX  (also,  see  IA  and  DAUX): 

TEST  =  -1.,  when  the  first  IB  equations  are  being  integrated 
(at  a  smaller  step-size)  in  the  multi-step  Cowell  mode. 

TEST  =  +1.,  when  the  N  equations  are  being  integrated  in  the 
Runge-Kutta  or  Cowell  mode. 

IDH  an  indicator  switch  to  the  user  after  each  integration  step. 

IDH  =  1,  if  the  step  size  has  not  changed 

IDH  =  2,  if  the  step  size  of  all  N  equations  has  been 

reduced  to  h  =  h/(2*IR)  where  h  is  the  current 
step  size  of  the  first  IB  equations. 

IDH  =  3,  if  the  step  size  of  all  N  equations  has  been  doubled. 

IDH  =  4,  if  the  step  size  of  the  final  N-IB  equations  has  been 

doubled  (possible  only  in  multi-step  mode). 
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NTRY 


JHH 


JHD 


VMIN 


VMAX 


T 


a  special  option  to  simulate  multiple  entries.  The  user  must 
set  NTR.Y  prior  to  using  the  calling  sequence. 

NTRY  =■  1  Setup  entry  (store  all  initial  conditions  first). 

NTRY  =  2  Normal  Runge-Kutta/Cowell  integration. 

NTRY  =  3  Integrate  in  Runge-Kutta  mode  exclusively. 

(Note:  NTRY  =  1  must  be  used  prior  to  the  other  two 
values.  See  Section  II-B  for  further  details.  ) 

an  option  to  control  the  halving  procedure. 

JHH  =  1  Reduce  step  size  and  return  to  Runge-Kutta  mode. 

JHH  =  3  Return  to  previous  step,  restore  all  conditions 

(i.  e.  ,  t,  y-^,  and  y'.')  at  that  step,  reduce  step  size 
and  return  to  Runge-Kutta  mode. 

(Note:  See  Section  II-C  for  exception  to  JHH  -  3  option.  ) 
(No  longer  used  in  DEQ4) 

All  doubling  is  performed  in  the  Cowell  mode  through 
the  accumulation  of  alternate  steps. 

the  location  of  one  cell  used  by  DEQ4  for  the  halving  and 
doubling  tests. 

VMIN  =  10"*  ER/H^  (computed  internally,  varies  with  H) 

Initially  (prior  to  NTRY  =  1)  the  user  must  set  VMIN  £  1. 
if  TEST  £  2.(i.  e.  ,  multi-step  mode).  This  initial  value 
allows  the  user  to  control  the  doubling  procedure  for  the 
final  N-IB  equations  if  the  multi-step  mode  is  used.  A  - 
larger  value  of  VMIN  (initially)  will  reduce  the  accuracy 
requirements  for  the  final  N-IB  equations  and  thus  allow 
them  to  be  integrated  at  a  larger  step-size.  A  value  of 
10.  £  VMIN  £  10-5  is  recommended. 

the  location  of  one  cell  used  by  DEQ4  for  the  halving  and 
doubling  tests.  VMAX  =  10-3  ER/H2  (  computed  internally, 
varies  with  H), 

the  location  of  the  independent  variable  tn  --  stored  initially 
by  the  user  and  incremented  automatically  by  the  subroutine 
during  the  integration  procedure.  T  may  be  reset  to  a  pre¬ 
vious  step  by  the  subroutine  during  the  starting  procedure 
or,  during  the  halving  procedure  when  JHH  =  3. 


-6- 


H  the  location  of  step-size  h  stored  initially  by  the  user 

and  modified  automatically  by  the  subroutine  during  the 
integration  procedure.  Initially,  DEQ4  sets  H  =  H/IR  as  the 
initial  step-size  for  both  Runge-Kutta  and  Cowell.  If  H  =  0, 
DEQ4  sets  H  =  0.  01  (H  can  be  positive  or  negative).  In  the 
multi-step  mode,  H  will  contain  the  current  step  size  of  the 
first  IB  equations. 

Y  the  location  of  N  dependent  variables  yi  --  stored  initially 

by  the  user  and  modified  automatically  by  the  subroutine 
during  the  integration  procedure. 

YP  the  location  of  N  first  derivatives  yl  --  stored  initially  by  the 

user  and  modified  automatically  by\he  subroutine  during  the 
integration  procedure. 

Y2P  the  location  of  N  second  derivatives  y'.1  --  computed  and  stored 

by  the  DAUX  subroutine  supplied  by  the  user.  If  the  two-group, 
multi-step  mode  is  used,  then  DAUX  should  be  programmed  to 
skip  the  evaluation  of  the  final  N-IB  second  derivatives  when 
TEST  is  negative. 


T 1  through 
T7,  F2P,  F  IP, 
and  DLT1 
through  DLT8 


are  temporary  storages  used  by  the  subroutine.  Each 
of  these  temporary  storages  must  have  a  dimension 
of  N  and  must  be  preserved  throughout  the  integration 
procedure. 


B.  MULTIPLE  ENTRIES 


For  purposes  of  compatibility  between  machines  and  different  versions 
of  Fortran  IV,  the  parameter  NTRY  is  used  to  simulate  multiple  entires. 

This  requires  the  user  to  change  the  value  of  NTRY  prior  to  using  the 
calling  sequence.  This  is  done  as  follows: 

NTRY  =  1,  Setup  Entry 

The  user  must  set  NTRY  =  1  and  store  all  initial  conditions  and  options 
prior  to  using  the  calling  sequence.  This  includes  N,  II,  12,  IB,  IR,  ER, 
HMIN,  HMAX,  YMIN,  TEST,  NTRY,  JHH,  VMIN,  T,  H,  Y,  and  YP. 

After  setting  NTRY  =  i,  the  user  must  CALL  DEQ4(N,  II,.  ..  ,DLT8) 
for  the  setup  entry.  The  subroutine  performs  various  tests  on  parameters 
(e.  g.  ,  sets  H  =  0.  01,  if  H  =  0),  sets  initial  switches  and  options  for  the 
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integration  procedure,  calls  DAUX  to  form  yl1  at  tg,  and  then  returns  to  the 
user.  No  integration  is  performed  during  the  setup  entry.  The  above  entry 
must  be  performed  prior  to  any  other  entry  (one  time  only). 

NTRY  =  2,  Integration  Entry  -  Runge-Kutta/Cowell  Mode 

This  is  the  normal  integration  entry  and  should  be  programmed  within 
a  DO  loop  (to  control  situations  such  as  maximum  steps  in  case  of  an  error). 
After  setting  NTRY  =  2,  the  user  must  CALL  DEQ4(N,  Ii, .  .  .  ,  DLT8)  to 
integrate  one  step  in  the  Runge-Kutta/Cowell  mode.  The  subroutine  inte¬ 
grates  one  step  (H/IR  initially)  and  returns  to  the  user.  The  first  eight 
steps  will  be  in  the  Runge-Kutta  mode;  thereafter  the  subroutine  will  inte¬ 
grate  in  the  Cowell  mode  (h  =  H/IR),  doubling  approximately  every  16  steps 
until  h  -*■  H.  If  the  two-group,  multi-step  mode  is  used  (TEST  5  2.initially) 
then  the  step  size  of  the  final  N-IB  equations  may  continue  to  double,  subject 
to  the  restrictions  of  the  initial  values  of  TEST  and  VMIN.  In  this  case, 
during  each  entry  to  DEQ4,  the  first  IB  equations  will  be  integrated  R  steps 
(h  =  H)  while  the  final  N-IB  equations  will  be  integrated  one  step  (h  =  R*H). 

All  of  this  is  automatic,  and  the  user  need  only  loop  through  the  calling 
sequence  to  integrate  from  t^  to  The  user  should  test  T  after  each 

integration  step  to  determine  the  proper  exit  procedure,  using  NTRY  =  3 
to  integrate  to  a  specific  time  point  If  output  is  desired  at  a  particular 

time  point  t^  (t.  ^  tn),  an  interpolation  subroutine  such  as  RW  NTRP  is 
recommended. 

NTRY  =  3,  Integrate  in  the  Runge-Kutta  Mode 

After  any  integration  step,  or  after  setup  entry,  the  user  can  integrate 
one  or  more  steps  in  the  Runge-Kutta  mode.  After  setting  NTRY  =  3  (and 
changing  the  step-size  H  if  necessary),  the  user  must  CALL  DEQ4(N, II,.  .  .  , 
DLT8)  to  integrate  one  step  in  the  Runge-Kutta  mode.  The  step  size  can  be 
changed  after  each  Runge-Kutta  integration  step,  and  negative  H  is 
permissible.  NTRY  =  3  should  be  used  primarily  to  end  the  integration  at 
a  specific  value  tj^ND'  ^  can  a^so  use<^  to  integrate  to  a  specific  output 
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time  fcouTPUT  "  H  <  fcOUTFUT  <  fcn  +  where  tn  is  the  current  value 
of  the  independent  variable  and  H  is  the  current  step  size,  if  the  regions 
TEST,  T,  H,  Y,  YF.  Y2P,  and  T1  through  T4  are  saved  and  restored 
before  continuing  the  normal  integration  (i.  e.  ,  NTRY  =  2).  Note  that  the 
step  size  for  the  Runge-Kutta  integration  should  be  <  H/1R  if  the  same 
accuracy  is  desired,  where  H  is  the  current  step  size  of  the  Cowell  mode. 

If  the  two-group  multi-step  mode  is  used,  then  several  Runge-Kutta  steps 

may  be  required  to  integrate  from  tn  to  fcouTPUT  e*  »  fcn  “  <  ^UTPUT 

<  t  +  R*H).  (See  Section  II-C  for  further  details.  ) 

C.  TWO-GROUP  MULTI-STEP  MOPE 

The  subroutine  assumes  that  the  N  equations  can  be  separated  into  two 
groups  (IB  and  N-J.B)  such  that  the  first  IB  equations  are  not  dependent  on 
the  final  N-IB  equations.  This  would  be  the  case,  for  example,  for  a 
system  of  nonlinear  differential  equations  which  are  integrated  simultaneously 
with  the  variational  equations  (i.  e.  ,  Newton-Raphson  Method).  Since  the 
variational  equations  are  always  linear  and  also  require  less  accuracy,  they 
can  be  integrated  at  a  larger  step  size  than  the  first  IB  nonlinear  differential 
equations.  This  is  accomplished  in  DEQ4  through  the  use  of  the  initial 
values  of  TEST  and  VMIN.  TEST  controls  the  maximum  ratio  of  the  two 
step  sizes  and  VMIN  controls  the  accuracy  of  the  final  N-IB  equations.  For 
example,  if  TEST  =  3-and  VMIN  =  1000.  initially,  then  the  maximum  possible 
ratio  of  the  two  groups  will  be  =  2^  =  2**(TEST  -1.)  and  the  accuracy 

requirements  for  the  final  N-IB  equations  will  be  ER*1000,  where  ER 
determines  the  accuracy  of  the  first  IB  equations. 

The  doubling  procedure  works  as  follows:  Initially,  all  N  equations 
are  integrated  in  the  Runge-Kutta  mode  at  a  step-size  H  =  H^/IR.  After 
eight  Runge-Kutta  steps  all  N  equations  are  integrated  in  the  Cowell  mode 
at  a  step-size  h  =  Hq/IR.  Thereafter  (assuming  halving  is  not  required),  the 
step  size  for  all  N  equations  will  double  every  16  steps  (in  the  Cowell  mode) 
until  the  step  size  for  the  first  IB  equations  can  no  longer  be  doubled.  From 
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that  point  the  step  size  of  the  final  N-IB  equations  will  continue  to  double 
under  the  restrictions  of  TEST  and  VMIN  described  above.  For  example, 
during  an  entry  to  DEQ4  (NTRY  =  2)  if  R  =  4  the  following  integration  will 
be  performed  in  the  Cowell  mode: 


FIRST  IB  EQUATIONS 
FINAL  N-IB  EQUATIONS 


H  [  H  |  H  |  H 
R*H 


fcn+l 


Note  that  during  one  entry  into  DEQ4,  the  first  IB  equations  are  integrated  R 
steps  with  a  step  size  of  H  while  the  final  N-IB  equations  are  integrated  one 
step  with  a  step  size  of  R*H.  Upon  the  next  entry,  if  halving  is  required  for 
the  first  IB  equations,  then  H  is  reduced  to  H  =  H/(2*IR)  for  all  N  equations 
and  the  above  procedure  is  repeated.  If  JHH  =  1  the  integration  will  continue 
from  T  =  tn+j.  If  JHH  =  3  the  integration  will  normally  continue  from  T  =  tn 
with  all  conditions  restored  at  that  point.  However,  if  the  final  N-IB  equations 
were  being  saved  for  doubling  at  T  =  tn  (i.  e.  ,  R^vIAX  had  not  been  reached), 
then  the  integration  will  be  continued  from  T  =  tn+j  (same  as  JHH  =  1). 

Thus,  the  doubling  of  the  final  N-IB  equations  is  given  precedent  over  the 
halving  procedure  in  regards  to  the  JHH  =  3  option.  For  this  reason 
TEST  £  5. is  recommended,  since  the  doubling  procedure  for  the  final  N-IB 
equations  is  discontinued  when  Rj^y^;  has  been  reached. 

Since  2(R-1)(N-IB)  derivatives  are  no  longer  computed  during  each 
Cowell  step,  the  use  of  the  two-group  multi-step  mode  will  result  in  a 
considerable  saving  of  computing  time  where  applicable. 

D.  CODING  INFORMATION 

The  calling  sequence  can  be  simplified  to 

CALL  DEQ4 

by  placing  all  parameters  of  the  calling  sequence  in  (labeled  or  blank) 
COMMON  and  using  the  name  DAUX  (identical  with  the  name  used  internally 
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by  DEQ4)  for  the  auxiliary  subroutine.  In  additon,  all  parameters  (i.  e.  ,  Y, 
YP,  Y2P,  etc.  )  must  be  dimensioned  within  subroutine  DEQ4* ,  and  the 
SUBROUTINE  statement  in  DEQ4  must  be  changed  to 

SUBROUTINE  DEQ4 

E.  SPACE  REQUIRED 

In  addition  to  the  parameters  in  the  calling  sequence,  approximately 
3571g  cells  are  required. 

Note:  The  coefficients  required  for  the  Cowell  integration  mode  are 
stored  in  DEQ4  as  20  digital  octal  constants  through  the  use  of  the  DATA 
declaration.  These  coefficients  were  formed  in  double  precision  on  the 
CDC  6600  and  are  the  octal  equivalent  of  the  60-bit  real  floating  point 
constants  of  the  CDC  6000  series  machines.  See  Appendix  A  for  a  com¬ 
plete  description  of  the  mathematical  method. 
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APPENDIX  A 


MATHEMATICAL  METHOD 


The  ASC  DEQ4  subroutine  [Floating  Point  Cowell  (Second-Sum) 
Runge-Kutta  Integration  of  Second-Order  Ordinary  Differential  Equations]  is 
prepared  to  solve  the  system  of  equations  defined  as 


y- '  -  ^(t.  y^ ....  yN.  yj. ....  yfj)  U  =  i  >  2, . . . ,  n) 
yi(to)=yio'  y'i^o^y'io  <i=  *>2 . N> 


(A-  1 ) 


In  case  none  of  the  L  involve  the  first  derivatives  yf,  time  is  saved  by 
indicating  this  in  the  set-up  (i.e.,  set  II  =  -1).  DEQ4  includes  a  fourth  order 
Runge-Kutta  subroutine  that  is  used  for  the  starting  procedure  and  for  the 
halving  procedure.  The  Runge-Kutta  subroutine  can  also  be  used  independently 
of  the  main  subroutine  through  the  use  of  NTRY  =  3.  For  the  sake  of  com¬ 
pleteness  the  Runge-Kutta  equations  are  given  in  the  following  paragraphs. 

A.  1.  RUNGE-KUTTA  EQUATIONS 

Let  y.  and  y'.  be  the  values  of  y.  and  y'.  at  t  =  t  :  f.  be  the  second 
m  m  '  i  1  jl  n  m 

derivative  of  y^  at  t  =  tn;  and  h  be  the  increment  (step  size)  of  the  independent 
variable  t.  The  Runge-Kutta  formulas*  used  in  this  subroutine  are  as 
follows: 


*J.  B.  Scarborough,  Numerical  Mathematical  Analysis,  Third  Edition, 
John  Hopkins  Press,  Baltimore,  Maryland  (1955)  pp.  300-301. 
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» 


k..  =  hf.(t  ,  y.  ,  y'.  ) 
ll  1'  n  1  in  ■'in' 


k., 


ki2  =  hfi(‘n  +T  yin+Iyi„+tkU'  yin  +  4) 

,  h  .  h  ,  ,  h .  ,  ,  i2\ 

ki3  =  hfi(tn+2’  yin+  2yin  +  8kii,yin+~rj 

I* 

ki4  =  hfi(tn  +  h’  yi„  +  hyL  +  2ki3’  yin  +  ki3> 

4yin  =  k[yin+i(kil  +ki2  +  ki3>] 

Ayin=l[kil  +  2ki2  +  2ki3+ki41 
yi,  n+1  “  yin  +  Ayin 
yU+i  =yin+Ayin 


(A-2) 


For  the  special  second-order  equation,  you  have 


y'/  =  R(t,  y^,  .  .  .  ,  yn)  (1st  derivatives  missing) 


(A-3) 


It  should  be  noted  that  k^  =  k.^,  so  that  above  formulas  reduce  to  the 
following  Runge-Kutta  formulas: 

k. ,  =  hf.(t  ,  y.  ) 
ll  1'  n 

ki2  =  h£i(ln  +  2’  yia  +  2yin  +  8kil) 
ki3  =  ui(‘n  +  h'  yin  +  hyi„  +  lki2) 

Ayi„  ”  h[yin  +  I  (kil  +2ki2>] 

(cont.  ) 
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i 


Ay?  =  4-  [k..  +  4k. 9  +  k.,] 

7in  6  il  i2  i3J 


yi,  n+1  ~  yin  T  Ayin 


y-  ,  <  =  y?  +  Ay? 
*1,  n+1  7m  *' 


in 


(A-4) 


where  k^  of  Eq.  (A- 2)  is  now  k.g  of  Eq.  (A-4). 

The  subroutine  can  be  made  to  take  advantage  of  this  fact  by  a  simple 
change  in  the  calling  sequence  (i.  e.  ,  set  II  =  -1)  and  thus  save  one  deriva¬ 
tive  per  integration  step  in  the  Runge-Kutta  mode. 

A.  2  COWELL  STARTING  PROCEDURE 

A  Cowell  (Second-Sum)  method  based  on  eight  differences  is  used  to 
continue  the  integration  after  the  first  eight  Runge-Kutta  steps.  The  user 
must  ask  for  each  integration  step  and  the  subroutine  will  follow  this 
sequence: 

a.  Perform  eight  Runge-Kutta  steps  of  size  h  =  Hstart/^ 

(one  Runge-Kutta  step  for  each  entry  to  DEQ4)  to  obtain 

yu,  y:r  yJJ  thiough  yig.  y:8,  yj^. 

b.  For  each  of  the  N  equations,  that  part  of  the  difference  table 
above  the  diagonal  line  is  constructed  in  an  accumulative 
manner  during  the  eight  Runge-Kutta  steps  (DEQ4  returns  to 
the  user  after  each  Runge-Kutta  step). 

To  start  the  Cowell  integration  we  need  "Fjjq,  '  F^g,  y^g,  y?g,  y?g, 

Ai7’  Ai6’  Ai5  ’  Ai4  ’  Ai3’  Ai2  '  Ail  ’  and  Ai0  *  The  values  yi8’  yi8’  and 
y?g  are  available  after  the  eight  Runge-Kutta  integration  step.  The  other 
values  are  computed  as  follows: 


'Fi9  a  yi4/H  +  W0y^0  +  WlyU  +  W2yi2  +  W3yi3 

+  W4y^  +  W5y-  +  W6yJ£  +  W?y''  +  WgvJJ 


(A- 5) 
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II 

o 

•H 

+  5yk/H  +  VJo  +  vi4I  +  Vh 

+  V3yi3  +  Vii  +  ViS  +  V6^6  +  V7y«  +  V84s 

(A-6) 

4  =  yi8  -  47 

(A- 7) 

4  =  4s-2yn  +  46 

(A- 8) 

4  =  y'it  - 3  4t  +  3  46  -  4s 

(A-9) 

4  =  4s  - 4  4? +  6  y'.6  - 4  4s +  44 

(A- 10) 

4  *  4s  -  Syi7  +  10yi6  -  l0yi5  +  544  -  yi3 

(A- 11) 

.VI 

A.~ 

i2 

=  48  -  647+1546  -  2045  +  1544  -  643  +4? 

(A- 12) 

a  VII 

Ail  =  ^8 

-  7y"  +  21  y:^  -  35ylg  +  35  y"  -  21  ylj  +  7  y "  -  y” 

(A- 13) 

4o“  =  4s  - 8  yn 1  28  46  -  56  45 

■+  70  y"  -  56y"  +  28  y"  -  8yl'  +  y" 

(A- 14) 

Before  going  lc  a  Cowell  step,  the  step-size  h  =  Hstart/IR  is  tested. 
Only  the  first  IB  equations  are  used  to  test,  where  1  £  IB  <  N.  For  the  first 
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IB  equations  we  determine 


V  = 


max 

i  <  i  <  IB 


max(|y.g|,  YMIN) 


(A- 15) 


3  2 

If  V  >  10  *  ER/h  ,  then  the  ratio  of  the  seventh  difference  to  function 

is  too  large.  Therefore,  h  is  reduced  to  h  =  h/(2#IR).  the  initial  conditions 
(i.  e.  ,  tQ,  y^Q,  y.Q  ,  and  y^ )  are  restored,  and  the  integration  procedure  is 
restarted.  This  procedure  is  repeated  until:  (a)  V  <  10^  *  ER/h^  or 
(b)  h  <  HMIN.  For  (a)  the  routine  proceeds  to  a  Cowell  integration  step 
with  the  current  value  of  h.  For  (b)  the  routine  sets  TEST  =  13.  and  returns 
to  the  user  --  indicating  a  programming  error  or  that  ER,  Hgj.  or  HMIN 
are  too  restrictive.  The  constant  YMIN  (equals  1.  if  unspecified)  prevents 
division  by  y  near  a  zero;  for  example,  in  the  siue  calculation  YMIN  =  0.  01 
avoids  difficulty  near  180  deg.  The  value  ER  (equals  IE  -  1 1  if  unspecified) 
allows  a  larger  h  if  chosen  larger,  say  ER  =  IE  -  8. 

A.  3  COWELL  INTEGRATION 
3  2 

If  Y  <  10  *  ER/h  after  the  eighth  Runge-Kutta  step,  we  begin  the 

Cowell  integration  with  predictions  of 


?  /  T  TT  TTT 

yi9  ’  hTFil0  +  N0>'i8  +  Nl4  +  N2Ai6  +  N3  Vs' 


.IV  ,  „  .V  ,  „  .VI  ,  kT  Avn 


+  N4A.^  +  N5A.3  +  N6Ai2  -r  n7^il 


N_  A. ,  +  N 


a  VIII  \ 
8AiO  ) 


(A- 16) 


yi9 


-■n  -•  h('Fio  +  :>0y^  +  N.A*,  +  N,aH  +  N.A?1 


i9 


1  17 


2i6 


3  15 


.  KT  A  IV  .j.  Av  .  •  A  VI  ,  JH  .  A  A  VIII  \  .. 

+  N4A.4  +  N5Ai3  +  N6Ai2  +  N7Au  +  NgA.0  )  (A- 17} 
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These  equations  use  the  row  of  the  difference  table  above  the  diagonal 

line;  only  this  row  is  needed  for  a  Cowell  step  and  is  kept  up  to  date  as  the 

integration  proceeds.  The  prediction  for  y£g  is  omitted  if  the  option  II  -  - 1 

is  used.  Now  from  y.g  and  y!g,  we  obtain  y!L  (from  DAUX)  and  then  com- 

17  17  VIII 

plete  the  row  of  differences  out  to  A..  under  the  diagonal  line  in  Table  A-l. 
I  XX  ^  x 

For  example:  A!g  =  y"  -  y(g  ;  =  A?g  -  At-,  ;  and  so  on. 

With  this  row,  we  calculate  the  corrected  values: 


%  *  h2("Fuo  +  'Vw  +  +  B2i“  +  b3a“ 


X  T*  AIV  ■  n  AV  X  H  AVI  X  n  AVI1  X  U  AVIII\  /  A  I  OX 

+  B4Ai5  T  B5Ai4  +  B6Ai3  +  B7Ai2  +BQAi{  )  (A- 18) 


y'n  *  h('Fw  +  +  b2a“  +  b3a“! 


•  IV  ’  V  VI  VII  •  VTTT\ 

+  B4AiS  +  B5Ai4  +  B6Ai3  +  B?4  +  B8*il  J 


(A- 19) 


From  these  we  get  corrected  values  for  yfg  (from  DAUX)  and  recalculate 
the  entire  row  under  the  diagonal  line.  Using  these  new  values  we  calculate 
new  values  for  y^  and  yf  using  Eqs.  (A-18)  and  (A-19).  Thus,  y.g  and  yfg  are 
corrected  twice;  however,  new  values  for  y'.'  are  not  computed  for  these  new 
values.  Next  we  compute  'F-^g  =  'F^g  +  y/'  and  /7F^  ^  =  ^F.jq  +  'F.jq.  This 
completes  the  integration  step  and  DEQ4  exits  to  the  user. 

On  the  next  entry  we  repeat  the  test  for  halving  [  i.  e.  ,  at  Eq.  (A-151j. 

If  halving  is  not  required,  we  continue  the  Cowell  integration  using  Eqs.  (A-16) 
through  (A-19)  as  before.  If  halving  is  required  after  at  least  one  Cowell 
integration  step  has  been  completed,  we  set  H  =  H/(2*IR)  and  start  a  new 
sequence  of  eight  Runge-Kutta  steps.  If  JHH  =  1,  we  do  not  retrace  ground; 
however,  if  JHH  =  3  (and  the  final  N-II3  equations  are  not  being  saved  for 
doubling  in  the  two-group  multi-step  mode)  then  we  return  to  the  previous 
step,  restore  the  conditions  at  that  point  and  then  start  a  new  sequence  of 
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eight  Runge-Kutta  steps  with  H  =  H/(2*IR).  This  can  be  shown  graphically 
as 


JHH  =  1 


current  step  (test  for  halving) 

H  /  set  H  =  H/(2  *  IR) 

,  COWELL  |  start  Runge-Kutta  step  from  t 

I _ I  i  i  i  i  i  t  i  i 

t  .  t  ^t  ,  4 

n-i  n  n+i 


JHH  =  3  [ 

t 


current  step  (test  for  halving) 
COWELL  \ 


n- 1 


n 


return  to  t  . ,  restore  conditions  (t  =  t  . ,  y . 
i  n-  x  ii ••  x  n-  x 


n-1  n 


^i,  n- 1  ’  ^i,  n-  1^ 


set  H  =  H/(2*IR) 

start  Runge-Kutta  step  from  tn  ^ 


During  the  halving  procedure  all  N  equations  return  to  Runge-Kutta 
with  the  same  step  size  (i.  e.  ,  the  two-group  multi-step  mode  is  temporarily 
suspended). 

If  halving  is  not  required,  we  proceed  with  a  Cowell  step,  making  the 
following  test  every  other  step  to  determine  whether  we  may  be  able  to 
double  h.  If 


v-  (10"2>7*er 

h2 

we  may  be  able  to  double  h. 

We  test  further  to  see  that 


(A- 20) 


W  = 


max 

1  <  i  <  IB 


VIH  1 

|  iO  |  <JJT_*_E R 

max  ( |  yig  | ,  YMIN)  ~  h2 


(A-Zl) 


Xf  the  conditions  of  Eqs.  1A-20)  and  (A-21)  are  satisfied  for  16  steps  (testing 
every  other  step)  then  the  step-size  h  is  doubled  for  all  N  equations  in  the 
Cowell  mode.  In  addition,  all  of  the  data  above  the  diagonal  line  is  doubled 
in  the  Cowell  mode  using  Eq.  (A- 23).  We  then  proceed  to  a  Cowell  step 
using  Eqs.  (A- 16)  through  (A- 19)  with  step-size  2H.  On  the  next  entry,  the 
halving  and  doubling  tests  are  continued  as  before. 

The  doubling  procedure  requires  a  counting  device  and  it  works  as 
follows.  Initially,  we  set  K  =  1  at  the  beginning  of  the  first  Cowell  step. 

The  counter  is  advanced  every  other  step  if  Eqs.  (A-20)  and  A-21)  are 
satisfied  for  the  first  IB  equations.  We  accumulate  the  data  required  for 
doubling  at  alternate  steps.  If  the  tests  fail  at  any  point  (assuming  halving 
is  not  required)  the  counter  is  restarted  if  all  N  equations  are  being  integrated 
at  the  same  step  size  (i.  e.  ,  TEST  =  VMIN  =  1.  initially).  However,  if  the 
two-group  mode  is  being  used,  then  the  doubling  tests  for  the  first  IB 
equations  are  temporarily  suspended  (i.  e.  ,  the  step  size  for  the  first  IB 
equations  remains  constant  unless  halving  is  later  required).  Thereafter, 
the  doubling  for  the  final  N-IB,  equations  is  continued  using  the  following 
test,  that  is,  if 


VAR  = 


max 


IB  <  i  <  N 


a  VH 

jfiLl 


[max  ( |  y.g  j  ,  YMIN)J 


10' 2  *  ER  *  VMIN 


(A- 22) 


for  the  remaining  steps  [K  =  1(1)9],  then  the  step  size  for  the  final  N-IB 
equations  is  doubled.  Also,  for  the  final  N-IB  equations  all  of  the  data 
above  the  diagonal  line  is  doubled  in  the  Cowell  mode.  Note  that  VMIN  of 
Eq.  (A-22)  is  the  initial  value  of  VMIN  for  NTRY  =  1.  The  counter  K  will 
be  set  to  one  for  the  next  entry,  and  the  doubling  procedure  will  be  restarted 
for  the  final  N-IB  equations. 
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The  equations  used  for  doubling  in  the  Cowell  mode  are  as  follows.  Let  A 
represent  a  difference  (i.  e.  ,  A*  =  A^,  A^  =  A?£)  at  step-size  h,  and  A 
represent  a  difference  at  step-size  2b.  Then  we  compute 


A1  =  2A1  -  A11 
A11  =  4A11  -  (4Am  -  AIV> 

Am  =  8Am-(12AIV-6AV+AVI) 

AIV  =  16AIV  -  (32AV  -  24AVI  +  8AVI1  -  A™) 

AV  =  5AIV  -  (10A111  -  10An  +  5A1)  +  (y"6  -  y£' ) 

AVI  =  6AV  -  (15AIV  -  20Am  +  15An  -  6a1)  -  (y"6  -  y" ) 

AVI1  =  7AVI  -  (21AV  -  35AIV  +  35A111  -  21  A11  +  7A1)  +  (y''6  -  y"  ) 

AVm  =  8AVI1  -  (28AVI  -  56AV  +  70AIV  -  56AUI  +  Z8AU  -  8A1)  -  (y"6  -  y" ) 
'F(A)  =  Yg/2h  +  WQ  Yq  +  WlY"  +  W2y"  +  W^' 

+  W4^8  +  W5*'i0  +  Wi  +  W7?14  +  *8*16 


#/F(A)  =y8/(2hr  +5yJ/2h  +  V0y^  +  y"  +  V2y" 

+  v3n' +  w +  v5yio + v6n'2  +  v7n4 + 


(A- 23) 


where  the  values  yg,  yg,  and  y^' ,  y2  , . .  •  ,  y'^'g  are  the  integrated  values 
saved  at  alternate  steps  of  the  Cowell  integration  procedure.  For  the  two- 
group  mode  h  will  be  the  value  in  cell  H  when  the  step  size  for  all  N  equations 
is  being  doubled  and  R  *  H  when  the  step  size  for  the  final  N-IB  equations 
is  being  doubled. 
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A.  4  CALCULATION  OF  COEFFICIENTS 


In  a  previous  version  of  the  subroutine  (ASC  DEQ2)  "F.jq  and  'F^  of 
Table  A- 1  were  calculated  by  first  forming  "F.j.  and  'F.,.  as  follows: 


'F._  =  y'.jH  -  Drt  y'.\  -  D.A*  -  D-A1*  -  D.A**1 
i5  yi4  07i4  1  i4  2  i3  3  i3 

n  aIV  n  aV  n  aVI  n  aVI1  n  aVI11 
"  D4Ai2  “  D5Ai2  “  D6Ail  “  D7Ail  "  D8AiO 


/it2  n  n  a11  n  aIV  n  aVI  n  aVI11 

F.r  =y../H  -  C„y..  -  C..A.-  -  C.A.  A  -  C,  A. .  -  CqA.a 
x5  y  iA  0^14  2  i3  4  i2  6  ll  8  iO 


Then  "F.jq  and  'F^g  v^re  calculated  through  the  sequence  as 


'FiK+i  =  'FiK +  y'ik  (K  =  5’ 6 . 8) 


(A- 24) 

(A-  25) 


(A- 26) 


'F  =  "F  +  'F 

i  K+l  iK  iK 


(K  =  5,  6 . 9) 


(A-  27) 


However,  by  substituting  ordinates  (y )  for  differences  (A.^)  (e.  g.  ,  Ai4  = 
y!^  -  yf^;  A^  =  y^  -  2  y +  y!^)  in  Eqs.  (A-24)  and  (A-25)  and  through  the 
use  of  Eqs.  (A-26)  and  (A-27)  'F.g  and  "F.jq  can  be  defined  by  Eqs,  (A-5) 
and  (A-6).  Tables  A- 2  and  A-3  give  the  value  of  and  W^.  Table  A-4  gives 
the  values  of  the  other  coefficients. 

The  coefficients  for  V.  and  W\  were  formed  in  double  precision  on  the 
CDC  6600  and  are  included  in  subroutine  DEQ4  as  20  digit  octal  constants. 

For  example,  using  Table  A-3 


V?  =  (2  -  5D5  -  5D6  +  35D7  +  40D8  -  C6  +  8C8)  (A- 28) 
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Table  A-2.  Coefficients  for  for  Bq.  (A~5) 
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APPENDIX  B 

FORTRAN  IV  SOURCE  LANGUAGE  FOR  DEQ4 
(CDC  6400-6600) 


B-  1 


oooooooooooo 


ASC  DEQ4  (8TH  ORDER)  R.K. /GAUSS  JACKSON  (4-18-67)  J.F.  HOLT 


OEQ4OOO 1 


SUBROUTINE  0EU4 <N, II » 12 , 1 A  »  IB , IR ,ER t HMIN *  UMAX , YMIN, DAUX *TEST , IDH,  OEQ4OOO2 

1  NTRY,JHH, JHD,VMIn,VMAX,  DEQ4OOO3 

2  T»H»Y»YP,Y2P»T1»T2»T3*T4«  DEQ4OOO4 

375«T6,TT,  -  DEQ40005 

4  F2P,FlP,0LTl,DLT2,DLT3,0LT4,DLT5,0LT6,DLT7,DLTR)  DEQ4OOO6 

DATA  A0»AltA2»A3»A4»A5»A6,A7,A8/0l7l452525252B2525252,Ol7l45252525DEQ40007 
*252525252 » 017 1450421 0421 0421 042 »01 7144631463146314631, 017144441671DEQ40008 
*441 67 1441 ,01 7 144272727272727272 ,01 7 1441 42 124345450046 ,01714402477  ODEQ40009 
*525655446,017137642352047260270/  DEQ4OOIO 

DATA  APO, API »AP2»AP3»AP4,AP5,AP6,AP7»AP8/01 71 74000000000000000 »0170EQ400 11 
*166525252525252525, 017166000000000000000, 017165447644764476447, 017DEQ40012 
*165216161616161616,017165031250377565231»017164674152142466774»01?DE040013 
*164557436021207661 ,017164456716206452330/  DEQ40014 

DATA  B0,Bl,82,B3,B4,B5,B6,B7,B8/0l7l45252525252525252to00000000000DEQ400l5 
*000000000 *060673567356735673567, 060673567356735673567,06070041 0313DEQ400 16 
*556630410, O6o7oi44167l44l67l441, 060702347642525026560 ,o60703132207DEQ400l7 
*010257763,060703607613672523307/  DEQ4OOI8 

DATA  BPO,BP1»8P2,RP3,BP4,BP5,BP6,8P7»8P8/017174000000000000000»06oDE040019 
*632525252525252525, 060642525252525252525, q6o651 1751 1751 1751 175,060DEQ40020 
*653 1463 14631 46 J 146, 06 q66o543327 060 160543, 06 0662134066 134 066I 34, 06oDEQ4002l 
*663153172725033257,060663753007252122331/  DEQ40022 

DATA  WO,Wl»w2,w3,w4,W5,V(6,w7,w8/0l7045506l00743i07612,o60700423525DEQ40023 
*521 066740, 017124634350133463213,060633425400145555331, 017174000000DEQ40024 
*000000000. 01720421651777] 51 11 2?, 017177663070575106313,017204007354DEQ40025 
*252256711,017177776456357607156/  DEQ40026 

DATA  VO»V1,V2,V3,V4,V5,V6,V7,V8/01706677240653607546?,060653305233DEQ40027 
*301450700, 017145726271235706017, 060612435756631615770, 017P14634266DE040028 
*567637354, U1722426640671 1557615, Ol72l56363635i5R16l35,ol72l4022662DEQ40029 
*701160620,017177770713144246212/  DEQ4003o 

the  calling  sequence  can  be  simplified  by  putting  all  parameters  DEQ40031 

EXCEPT  NTRY  and  I.UUX  in  LABELED  common,  (SEE  EXAMPLE  BELOW)  DEQ40032 

HOWEVER*  THE  DIMENSION  STATEMENTS  MUST  BE  CHANGED,  OR  INCLUDED  DEQ40033 
IN  THE  LABELLED  COMMON  STATEMENT,  THE  NEw  CALLING  SEQ.  BECOMES  DEQ40034 
CALL  DE«4  <NrKY,DAUX)  DEQ40035 

ALSO,  THE  SUBROUTINE  STATEMENT  (SEc!  ABOVE)  MUST  BE  CHANGED  TO  DEQ40036 
SUBROUTINE  0Eu4 (NTRY, DAUX)  DEQ40037 

EXAMPLE  FOLLOWING  SHOWS  LABELED  COMMON  FOR  SIMPLIFIED  CALLING  SEQ.DEQ40038 
COMMON/CDEO/N»Il,T2»lA»IB»lR»ER»HMlN*HMAX,YMlN»TEST»IpH*IRK»JHH»  DEQ40039 

1  Jh0,VmIN,VMAX,T,H,Y,YP»Y2p,T1,T2,T3,T4,T5,t6,T7,  DEQ40040 

2  F2P,FlP,DLTl,0LT2,DLT3»  DEQ40041 

3  DLT4,DLT5,DLT6,ULT7,0LT8  DEq40042 
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C  DIMENSION  Y(  63)  ♦  YP  (  63)»Y2P<  63), Tl(  63),T2(  63>»T3(  A  3)  ? 

C  1  T4(  63) ,T5  (  63) » T6  {  63),T7(  63), 

C  2  F2P(  63) »F1P  (  63) »DLT1 (  63>*DLT2(  63), 

C  3  DLT3  <  63),DLT4(  63),DLT5(  63),dlt6<  63),DLT7(  63), 

C  4  DLTS<  63) 

Y(2) ,YP(2) »Y?P(2) »Tl (2) ,T2(2) »T3(2) ,T4(2) 

T5(2) ,T6(2) ,T?(2) 

P2P (2) »FiP(2> 

DLT 1 <2) ,DLT2(2) »OLT3(2) ,DLT4(2) ,DLT5(2) ,DLT6(2) 
l)LT7(2)  ,DLT8<2) 

OPTIONS  AND  INITIAL  CONDITIONS  (T,H,y»YP)  MUST  BE  STORED 
TO  THE  SETUP  ENTRY  TO  DEG4.  (DEQ4  ENTRY  CALLS  DAUX) 
SPFCIAL  NOTE  *«******#**********#*e***<H*#**#*#*  SPECIAL  NOTE 
WITH  THE  DATA  ALLOCATION  SCHEME  OF  0E04,  -IT  ?S  POSSIBLE 
THE  USER  TO  INTEGRATE  THE  FIRST 
COURSE,  THIS  ASSUMES  THAI  THE 

THE  FINAL  N-P  fQUATIONS.  (SEE  WRITE-UP  FOR  FURTHER 
IS  THE  NUMBER  OF  EQUATIONS. 

AN  OPTION,  SUCH  THAT  IF 

G.E.  +0,  Y2P  =  F<X,Y,YP)  U.E.?  FUNCTION  OF-  1ST 
L#T  «  0,  Y2P  =  F (X, Y) 

OPTION,  SUCH  THAT  IF 

G.E,  ♦  !)»  VARIABLE  STEP-SIZE  MODE  IS  USED. 

L.T,  n,  rlXEO  STEP-SIZE  MODE  IS  USED. 

INDICATOR  SWITCH  TO  THE  USER  DURING  EXIT  TO  DAUX  SUB. 
FOR  1ST  (1,2,3  IF  R.K.)  PASS  THRU  DAUX, 

FOR  FINAL  PASS  THRU  DAUX.  (APPLIES  FOR  EACH  STEP.) 

IN  THE  COWELL  MODE,  lA="l  WHEN  THE  DERIVATIVES  OF 
PREDICTOR  ARE  BEING  ASKED  FOR,  AND  lA=+l  WHFN  THE 
DERIVATIVES  OF  THE  CORRECTOR  ARE  BEING  ASKED  FOR. 

ONLY  THE  1ST  IB  (L.E.  N)  EONS.  ARE  TESTED  DURING  E^ROR  _ 

FOR  THE  TWO  GROUP  MODE*  THE  FIRST  IB  EONS  ARE  INTEGRATED  AT 
STEP-SIZE  H,  AND  THE  FINAL  N-IB  EQNS  ARE  INTEGRATED  AT  STEP- 
SIZE  R*H  (R=  1.,  2.1  4.,  8.,  ETC.  DEPENDING  ON  THF  INITIAL 
VALUES  OF  TEST  AND  VMIN).  (SEE  TEST  AND  VMIN) 

FOR  A  GIVEN  STEP-SIZE  H,  COWELL  STEP  =  H,  R.K,  STEP  =  H/IR. 
RELATIVE  ERRoR  CPITERIA. (SET  ER=l,E-9  OR  LESS) 

IS  THE  MINIMUM  STEP-SIZE  ALLOWED.  (ABSOLIJTF  VALUE) 

THE  MAXIMUM  STEP-SIZE  ALLOWED.  (ABSOLUTE  VALUE) 

THE  MIN,  aBS.  VALUE  OF  Y { I )  ALLOWED  FOR  THE  ERROR  TEST. 

IS  an  external  sub.  supplied  BY  USER  TO  eval.  the  2ND  deriv. 
AFTER  EACH  integration  step, 


dimension 
dimension 
dimension 
dimension 
dimension 
c*##**note  all 
prior 

c*#*» 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


of 

ON 

N 

II 


12 


IA 


IS 

II 

11 
IS 

12 
12 
IS 


IS 

IS 

AN 

IS 

IS 

AN 


IA=-1 
I  A=.  1 
NOTE 
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DEQ40043 
DEQ40044 
DEQ40045 
DEQ4Q046 
DEQ40047 
DEQ40048 
DEQ40049 
DEQ4OOSO 
DEQ40051 
DEQ40052 
DEQ40053 
DEQ40054 
#»##CDEQ40055 
FOR  DEQ40056 

P  EQUATIONS  SIMPLY  BY  CHANGING  N,DEQ40057 
FIRST  P  EQUATIONS  DO  NOT  DEPEND  DEQ40058 

DETAILS*)  DEQ40059 
DEQ4OO6O 
DEQ4OOGI 
DERIV* ) DEQ4OO62 
DEQ40063 
DEQ40064 
DEQ40065 
DEQ4OO66 
DEQ4OO67 
DEQ4OO68 
DEQ4OO69 
DEQ40070 
DEQ40071 
DEQ40072 

tests, DEQ40073 


THE 


IR 

ER 

HMIN 

HMAX 

YMIN 

DAUX 

TEST 


IS 

IS 


DEQ40074 

DEQ40075 

DEQ40076 

DEQ40077 

DEQ40078 

DEQ40079 

DEQ4OO8O 

DEQ4OO8I 

DEQ4OO82 

DEG40083 

DEQ4OO84 
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OOOOOOOCJUOOpOO 


* 


TEST  WILL  BE  *1.  IF  THE  INTEGRATION  WAS  A  RUNGE-KllTTA  STEP*  DEQ4OO85 
TEST  WILL  BE  -I.  IF  THE  INTEGRATION  WAS  A  COWELL  STEP,  DEQ4OO86 
NOTE  INITIALLY  SET  TEST*1,  FOR  NORMAL  INTEGRATION.  FOR  TWO-DEQ4OO67 


DEQ4OO88 
DEQ4OO89 
— OEQ40090- 
DEQ40091 
DEQ40092 
DEQ40093 
-  -  DEQ40094 


GROUP  MODE  SET  TEST  G.E.  2.  <MAX  R  *-r2,**TTEST-l, ) 

■  IDH^.AFTER  EACH  INTEGRATION  'STEP t  ‘ 

BE  >I-T^:TttE=1^^'-SiZ£--HAS--NO-f^cHANGEbv  -i  1  -  - 
I OH  WILL  BE  *2  IF  THE  STEP-SIZE  HAS  BEEN  HALVED. 

IDH  WILL  BE  *3  IF  THE  STEP-SIZE  HAS  BEEN  DOUBLED 

^-rtOH  WILL  BE  +4  -IF  H-EOR -FINAL  N-IB  EQNS.,  DOUBLED-.  -  - 

:-rV%  -NYRYL  is."*  SPECIAL  0  P  T I ON^T-Q :,ALL0 HU  LT IPLELE  N  t RiESY"^ 

i;^^MCtiE'Tf.;pT.8fe  -'^fliRi'fXNJUiALaeoND^ttOfe^a  entry.*.!.  .  j>e<$oo«&. 

- :-=:,=^=NY-RY:^*--N0RMAL^.>f.y&AU'SsVjACKS^N  INTEGRATION.  --  DEQ40096' 

NTRY=3*  INTEGRATE  IN  R.K.  MODE  EXCLUSIVELY.  DEQ40097 

JHH  IS  AN  OPTION  TO  CONTROL  THE  HALVING  PROCEDURE.  IF*  DEQ40098 

t;  -.  *.:■  -JHH*Ir-  a«IR  R*K.rS?f»S-ARET-LJSE&  DURING  THE- HALVING  PROCEDURE. DEQ4O 099 
•eg.  <  JHH»3»  return  Tot f  (Tt*  -"iLT TEH£T::Jo.'-JTI:+  I) '»  8*IR  RKDEQ40100. 

CL  ^7';^  [j_=SE£  HJLITE-UP  FOR.  PiiRTWERLog^I LS:  _£ON CERNJNG  Jfffi. LV.TN &  pR.nCEOUREDEQ4.Ol 0 1 
&-  =‘-OhD'  -IS  '^GPT-ION  W-OO'NTROL^fHE^OOUBL-ING  'PRbcEOLlRE.  ^ *  - ' "  *  DEQ4OIO2- 

c  NOTE  JHD  IS  NO  LONGER  Used  in  DEQ4.  all  DOUBLING  is  done  DEQ40103 

c  IN  THE  COWELL  MODE  BY  THE  ACCUM.  OF  ALTERNATE  DATA.  OEQ4OIO4 

C;-_-Z->^HIN  IS  THE--  LO.C  OF^  l;  CELL  .  U^Or  &Y  £DEQ.  FOR  HALVING/OOUbUNG  TESTDEQ4OIO5 
ROt^  .  iNl T TALLY  SET  VM in* lV£P  fP«2.»3,»,..5);:FOR  TWO-GROUP  HODEDEQ4OIOE 
Cssgr^VNA X  YSLEHErlLOC  OF  J.  CELL  USED  BY  pDEQ  FOR.  HAlVING/doUbLING  TESTD.EQ4OIO7. 

c*- ™^^TS-^’E:-rocATiON  or  the -independent  VARiABLrTar.  ~7~:  DEoioiog 

C  H  IS  THE  LOCATION  OF  THE  INITIAL  AND  CURRENT  STEP-SIZE.  DEQ4OIO9 

C  Y  IS  THE  LOCATION  OF  N  DEPENDENT  VARIABLES.  DEQ4OIIO 

C-  .  -YP  -  IS  THE.  LOCATION  OF  N  1ST  DERIVATIVES.  :  DEQ4OIH 

C  Ygp  -  IS  THE  LOCATION.  OF.  N  2ND  .DERIVATIVES.  (COMPUTED  BY  DAUX)  DEQ4OH2 

THRU  T.7  ARE  VECIOR&  OF::Ol’MENSlQN  N  USED  BY  THE  SUBROUTINE.  DEQ4OII3 

Ci  =— f 2P.  ‘FT Py  A ND  DLT-l^THRU-Dt'TB'ARF 'VECTORS'  OF-  DIMENSION  N.  DEQ4OIU 

C  SETUP  FOR  DEQ4  SUBROUTINE.  CALL  DEQ4 (N, H , 12, . . . , DLT7*DLT8)  DEQ40115 

GO  TO  (1*2*3)  »NTRY  DEQ4OH6 

2  -  IDHbI  -  c  v  -----  --  -  DEGS4OI 17 . 

RXH*TEST,  *L.  t  OEQ4OII8 

-IF(RXH-i. >380*331. 381-  V-*  :  _  DEQ40T19, 

380  "~RXH»1 .  *  .  DEQ40'r20 

381  IKR*1  DEQ40121 

IF(VMIN-1. >382*383*383  DEQ40122 

382  VMIN=1000.  ---Cu.*"  -  -  -  -*  ,  r  .  .  -  DEQ40123 

383  VARX*VMIN  Lj  .  .  DEQ40124 

.  .  RXAsRXH...  i,  ^  *.%•:  S .  *  -  :DEQ40125 

--s^s-:fesT«W*-™-  v  r>- — --=* -J  DEQ40l26 
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TST— ♦ 1  , 

NSTR'i  =-l 

IF(N)391, 390*392 

390  N=1 

391  N  —  I ABS (N) 

392  IF(N-TB) 393.394*304 

393  IB=N 

394  IF  ( IR) 400 ♦ 400 , 41 0 
400  IB*N 

410  IF  ( I W) 420 ♦ 420 ♦ 430 
420  I Rs 1 6 

430  IF(H)432, 431*432 

431  H= ,  0 1 

C431  H=.01D0 

432  R=IR 
H=H/R 
HA=H 

IF(ER)433*433,434 

433  ER= 1 . E- 1 1 
C433  ER=1 *0-11 

434  HH=H*H 
N1=N 

IF  (YMINJ439, 440*450 
440  YMIM=1. 

C440  YMIN=1.D0 
439  YMIN-  ABS(YMIN) 

C439  YMTN=0AHS(YMIn) 

450  IF  (HMAX)4b0*460*470 
460  HMAX=1. 

C460  HMAX=1 .00 

470  IP(HMIN)471.47l,472 

471  HMIN*l.E-5 
C471  HMIN=1.D~5 

472  VMIN=ER/ (HH*100.) 

VMAX=VMlN»100000. 

I  A=- 1 

C  IA=-1  FOR  1ST  (1,2,3  IF  R.K.jPASS  THRU  QAUX 

NHH=+1 
CALL  OAUX 

C  SAVE  INITIAL  CONDITIONS  FOR  RFSTART 
TZ=T 

DO  473  1=1, Ml 


0EQ40127 

DEO40128 

DEO40129 

DEQ40130 

DEQ40131 

DEQ40132 

DEQ4ul33 

DEQ40134 

DEQ40135 

DEQ40136 

DEQ40137 

DEQ40138 

DEO4O 139 

DEQ40140 

DEQ40141 

DEQ40142 

DEC40143 

DEQ40144 

DEQ40145 

DEQ40146 

DE040147 

DEQ40148 

OEQ4 0149 

DE040150 

DEQ40151 

DEQ40152 

DEQ40153 

DEQ40154 

DE040155 

DEQ40156 

OEQ40157 

DEQ40158 

DEO40159 

DEQ40160 

DEO4OI6I 

DEQ40162 

DEQ4OI63 

DE040164 

DE040165 

DEQ4OI66 

OEQ4OI67 

OEQ4OI68 
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f  1 
‘  - 

* 

T5 ( I ) »Y ( I ) 

DEQ4OI69 

T6 ( I ) =YP  ( I ) 

DEQA0170 

473 

T7 ( I ) =Y2P ( I ) 

DEQ40171 

' 

IF (HMAX-HMIN) 261 .474.261 

DEQ40172 

• 

474 

NHH*-1 

DEQ40173 

r 

V 

EXIT  FROM  DEQ2  SETUP 

DEQ40174 

GO  TO  261 

OEQ40175 

c 

normal  integration  entry.  CALL  RKCw 

DEQ40176 

c 

ENTRY  RKCW 

DEQ40177 

c 

INTEGRATE  1  STEP  IH  OR  H/IR)  IN  THE  COWELL/R.K. 

MODE 

DEQ40178 

\  ? 

N1»N 

OEQ40179 

: 

IF(TEST)900.500»500 

DEQ40180 

C 

R.K.  STEP  (bR  R.K.  STEPS) 

DEO4OI8I 

bOO 

IDH=1 

DEQ40182 

bOl 

IF(IKR)505»502,502 

DEQ40183 

c 

R.K.  STEP  (INITIAL  STARTING  PROCEDURE) 

DEQ4OI84 

1  b02 

K=1 

DEQ4OI85 

1 

TEST =T5T 

DEQ4OI86 

IKRa-1 

DEQ4OI87 

"  505 

GO  TO  (5l0»520»530f540»550»560,570»580.590) .K 

DEQ4OI88 

[i  C 

K*1 

DEQ40189 

i  510 

DO  511  1=1. Nl 

DEQ40190 

DLTI ( I ) SY?P ( I ) 

DEQ40191 

I 

F1P(I)*W0*Y2P(I) 

DEQ40192 

5X1 

F2P ( I ) =V0*Y2P ( I ) 

DEQ40193 

„ 

512 

K=K+ 1 

DEQ40194 

513 

IF(K-10)3. 600.3 

DEQ40195 

- 

C 

SPECIAL  R.K.  INTEGRATION  ENTRY,  (NTRY=3) 

DEQ40196 

C 

ENTRY  RK 

DEQ40197 

C 

INTEGRATE  1  STEP  (H,  OR  H/IR)  IN  THE  RUNGE-KUTTA  MODE* 

DEQ4O 198 

C 

NTRY=3  IS  A  SPECIAL  USAGE*  INTEGRATES  ONLY  IN 

R.K.  MODE 

DEQ40199 

1  3 

N1=N 

DEQ40200 

| 

H1  =  H 

DE040201 

| 

HH=H*H 

DEQ40202 

H2=H/2, 

DEQ40203 

H3=HH/2. 

DEQ40204 

H4=H3/2. 

DEQ40205 

H5=H4/2. 

DEQ40206 

H6=H/6. 

DEQ40207 

DO  100  1=1 ♦ Nl 

DEQ40208 

TllI)sY(I) 

OEQ40209 

■ 

’ 

, 

T2 ( I ) =YP ( I ) 

DEQ40210 

- 

i 

f 

l 

'\ 

B-6 

- 

- - - - 

T3 ( I ) =Y2P ( I ) 

100  T4 ( I ) *Y2P ( I ) 

110  T=T*H2 
IA=-1 

IF  (II) 300*200*200 
C  R.K.  STEP  FUR  Y?P=F(X»Y»YP> 

200  DO  210  1=1 »Nl 
C  Y2P*F ( X ♦ Y » YP)  IF  Il=+ 

Y ( I ) =T1 ( I ) *H2ttT2 ( T ) +H5#Y2P ( I ) 

210  YP(I)sT2(I)+H2*Y2P(I) 

CALL  OAUX 

DO  220  I=1*N1 

T3 ( I ) =T3 ( I ) +Y2P ( I ) 

T4!I)*T4(I)+2.«YHp(I) 

220  YP ( I ) =T2 ( I ) *H2*Y2P ( I ) 

CALL  DAUX 

DO  230  I  =  1 » N 1 

T3 ( I ) =T3 ( I ) + Y2P ( I ) 

T4(I)=T4(I) ♦2,*Y2P(I) 
Y(I)=T1(I)+H1*T2(T)+H3*Y2P(I) 

230  YP(I)=T2(I)+H1*Y2p(I) 

T=T+H2 
CALL  OAUX 
DO  240  IsitNl 
240  T4(I)sT4(I)*Y2P(I) 

250  DO  260  Ial.Nl 

Y(I)  =T1(I)+H1#(T?(I) +H6ttT3 ( I ) ) 

260  YP ( I ) =T2 ( I ) ♦Hb* ( T4 ( I ) ) 

I  A=>  1 

C  IA=+1  FOR  FINAL  PASS  THRU  OAUX 

CALL  DAUX 
262  TEST=TST 

261  RETURN 

C  R.K.  STEP  FUR  Y?P=F (X*Y) 

300  DO  310  I = 1 i N 1 
C  Y2P=F { X * Y >  IF  1 1=- 
310  Y(I)=T1 (I) ♦H2«T2(T) ♦H5ttY2P(I) 

CALL  DAUX 
DO  320  Ial.Nl 
T3 ( I ) =T3 ( I > ♦2.*Y?P(I) 

320  T4(I)=T4(I)+4.«Y2p(I) 

T=T+H? 


DE040211 

DEQ40212 

DE040213 

DEQ40214 

DE040215 

DEQ40216 

DEQ40217 

DEQ40218 

DEQ40219 

DEQ40220 

DEQ40221 

DEQ40222 

DEQ40223 

DEQ40224 

DEQ40225 

DEQ40226 

DEQ40227 

DEQ40228 

DEQ40229 

OEQ40230 

OEQ40231 

DEQ40232 

DEQ40233 

OEQ40234 

DEQ40235 

DE040236 

DEQ40237 

DEQ40238 

DE040239 

DEQ40240 

DEQ40241 

DEQ40242 

DE040243 

OEQ4024-4 

DE040245 

DE040246 

DEQ40247 

DEQ40248 

DEQ40249 

DEQ40250 

DEQ40251 

DEQ40252 
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DO  330  I*1»N1 

330  Y(I)*T1  (I)  ♦MIMS  (I)  ♦H3*Y2P(I) 

CAUL-  DAUX  .  f  ' 

'  -=z'.  ipO  340  I-1»N1  '  j_‘  \  •  "  '  j 

^M-'-T4iX)^r4:(:X^Y2PtI)  ^ I".-  ^  ■ 

ei=H=^%'0^-T;0~250‘“  ‘  "  ^  - ^ 

C  K*2 

520  DO  521  1*1 »N1 

^^4P1^I*«DLT  1 IX >~8<*Y2P:< ;••>  H  •v"- 

^PfS^^TT)Wf*YZRTT  - “ - 

60  TO  512 


C  K*3 

g530?-?^)0:-53-l_  I«l  »NI  -  -  •  * ;  .■ --‘.w:  *3?sas==7-.--r-  -i^.-  ’-^r— r-_  -- - 

( imn  a )  +.28.^Y2Pa )  753^  -  ^Z?  ;  -~^cSi 

Y-2Pit> -*m.xz±M  - -i  .cm :  ::  . 

■===i“^t^tT)*Y2Ptl:>-  -«-=-;— •*•-  -  * - s - 

F1P(I)*F1P(I}*W2*Y2P<I> 

531  F2P ( I ) *F2P ( I ) +V2*Y2P ( I ) 

]ja~  .  ==£0HTQ  Si2::  -  --? -j£~  rft  --h^l -'. - 

csripfr**'"  ?  r^T  I  .  <”'*  •'-  •  :-:  -•-•"f' 

^16^Pr^4Xi6rt  jN1~ ---  -  -•  -  '-  _i~  -v. 

^^tTlWOLTl  ( I )  -56 **Y2P (I)  ' ‘ ’  '  '  "  '  ”  ' 

DLT2 ( I ) *DLT2 ( I ) -21 .*Y2P { I ) 

DLT3 ( I >  *DLT3 ( I ) -6.*Y2P < I ) 

^X^0LT4 ( l)sY2P ( I )  " 

>”  .  i  fFlP(I>«FlP(I)*W3*Y2Pd).-.  -  1- 

g5»-lte2P4d'*F2P  ( 1  )-»va»;Y,2P-j^  ll  -  .  -:>-  -  ->  - 

'aswasvQQisp^  tffc-  -  •■■■'"•  -  -  *  -------  -  "•  ---  - 


c 

550 


551 

C-\  . 

=-560" 


K*5 

DO  551  I*1*N1 

DLT1 < 1 ) *DLT1 ( I ) ♦70.*Y2PXIi 
DLT2(I)*DLT2(I)*35.*Y2P(lr  ' 

:DU3  <  I ) *DL.T3  ( X )  ♦  15.*Y2P  (j)  ,  ,  ..-Xl-Wfci 

T)LT4<I)*5.*Y2P<I>..DIT4<I)  ‘  ----- 

DLT5(I)*Y2P(I) 

F1P ( I ) *F1P ( I ) +W4*Y2P ( I ) +YP ( I ) /HA 

F2P ( I ) *F2P (I ) ♦  V4*Y2R.U >* (&**YP (  j) /HA) MY t‘I ) / <HA*HA ) ) 


GO  TO  512 
K*fc 

;D0'-'56 1~-F*TVN  T 


DEQ40253 

DEQ40254 

-  OEQ40255-- 
. DEQ40256- 

;2.  i  i:?OEQ4.02-57| 
~  -  -=^-DEQ4025B= 
DEQ40259 
DEQ40260 
.  i  r  :^%OEQ4Q26_ll 

-  Vx---g:;-^?.EQ4t)2»2r 

—  - - - DEQV0264 

DEQ40265 
DEQ40266 
;--^0E«4-0267s 
-t  r:$DEQ40‘2685 
JlJ?  _  ^QEG4:02&9} 
.->.'^«s3^Eq^2YqF 
DEQ40271 
DEQ40272 
v  ™.nEQ4.0273i 

■"--S  ^  saEQ4027si 
^--qe^o^T^ 

DEQ40277 
DE040278 
-  ..-DEQ40279; 

-  DEQ4O280 

._  '-  -SOEG4Q281: 

'—  "-DEQ40282- 
DEQ40283 
DEQ40284 
0EQ40285: 
DEQ4O286; 
tX-=DEQ^(l28X 
DEQ40288- 
DEQ40289 
DEQ40290 
-  ^  DEQ40291 

v;  DEQ40292 

-  3DEQ4:029.3 

- 1  'DEQ40294 


DLT1<I)=DIT1 (I)-56.*Y2P(I) 

DEQ40295 

DLT2 ( I ) ~DLT2 ( I ) -35«#Y2P ( I ) 

DEQ40296 

DlT3<I)=DLT.3(I)-2n.*Y2P(I) 

DEQ40297 

DLT4 < I ) =DLT4 ( I ) -1 o .*Y2P ( I ) 

DEQ40298 

DLT5 ( I ) =DLT5 ( X ) -4 ,*Y2P < I ) 

DEQ40299 

DLT6 ( I ) =Y?P ( I ) 

DEQ40300 

F1P(I)*F1P(I)*w5*y2P(I) 

DEQ40301 

561 

F2P(I)=F2P<I) ♦V5*Y2P(I) 

DEQ40302 

GO  TO  512 

DEQ40303 

K=7 

DEQ40304 

570 

DO  571  1  =  1  *  N 1 

DEQ40305 

DLTKI)=DLT1  (I)+2fl.*Y2P(I) 

DEQ40306 

DLT2(I)=DLT2(I) +2l.*Y2P(I) 

DEQ40307 

DLT3 ( I >  =DLT3  < 1 ) ♦15.*Y2P(I) 

DEO40308 

DLT4  ( I )  =DI_T4  (I )  + 1  o  .*Y2P  ( I ) 

DEQ40309 

DLT5 ( I) =DLT5 ( I) +6 . *Y2P ( I ) 

DEQ40310 

D1.T6  ( I )  =3.*Y2P  ( I )  -DLT6  ( I ) 

DEQ4031 1 

DLT7(I)=Y2P(I) 

DEQ40312 

FlP(I)sFlP(I)+w6*Y2P(I) 

DEQ40313 

571 

F2P(I)=F2P(I)+V6*Y2P(I) 

DEQ40314 

GO  TO  512 

DEQ40315 

K=8 

DEQ40316 

580 

DO  581  1  =  1  *  N 1 

DEQ4O3I 7 

DLT 1 ( I ) -DLT 1 ( I ) -8«*Y2p ( I ) 

DEQ40318 

DLT2(I)=DLT2(I)-7f*Y2P(I) 

DEQ40319 

DLT3(I)=DLT3(I)-6.«Y2p(I) 

OEQ40320 

DLT4 ( I ) =DLT4 ( I ) -5,*Y2P ( I ) 

DEQ40321 

DLT5 ( I ) =DLT5 ( 1 ) -4 ,*Y2P { I ) 

DEQ40322 

DLTo(I)=[)LT6(1)-?.*Y2P(I) 

DEQ40323 

DLT7 ( 1 ) SQLT7 ( 1 ) -2.*Y2P ( I ) 

DEQ40324 

DLT8(I)=Y2P(I) 

DEQ40325 

F1P<I)*F1P(I) +W7*Y2P ( i ) 

OEQ40326 

581 

F2P ( I ) =F2P ( I ) +V7*Y2P ( I ) 

DEO40327 

GO  TO  512 

DEQ40328 

K=9 

DEQ40329 

590 

DO  591  I=l»Nl 

DEQ40330 

DLT1(I)=DLT1(I)+Y2P(I) 

DEQ40331 

DLT2  ( I )  =D1_T2  ( I )  *YpP(I) 

DEQ40332 

DLT3 ( I ) =DLT3 ( I ) ♦ Y?P ( I ) 

DE.Q40333 

DLT4(I)sDLT4(I)*y?P(I) 

DEQ40334 

DLT5 ( I ) =DLT5 ( I ) *Y?P ( I ) 

DEQ40335 

DLT6 ( I ) =DL 16 ( l ) +Y2P ( I ) 

DEQ40336 
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oo 


DLT7 (I) «DLT7 ( I ) +Y2P  ( I ) 
DLT8 ( I ) *Y2P ( I ) -DLT6 ( I ) 
F1P ( I ) *F1P ( I ) *W8*Y2P ( I) 
591  F2P(I)*F2P(I)+V8*Y2P(I) 

60  70  512 
--£NTRY  Cw:  - 

600  K=1 

KC»1 
■  -  KS*1 
RAT*i •' 

-v  £r:Mus»*V.~;  •  .  . . 

“  7'NR2”NHH‘  ~  :  ‘ 


NVR=1 
NJH3**1 

TEST  FOR-  HALVING  AND  DOUBLING- ;  . 

IDH=1  IF  STEP-SIZE  HAS  NOT  CHANGED •  iL. 

0  IDH-1  .  .  _  -f^Y; 

’  ICbIB+1  :  -------  ------- . --  - 

IF  (12)1000.901.901 
901  IF(NH2)9011» 9010. 9010 

9010  DO  902  1=1, IB 

IF (  ABS(0LT2(I)/AMAX1(  ABS <Y (I) ) ,YMIN) ) .GE.VMAX)  GO  TO  1100 
C  '  IF (DABS <0LT2 (I) /DMAX1 (DABS ( Y ( I) ) , YHIN) ) .G&iV.MAX)  GO. TO  1100 
■902  CONTINUE 

C  IF  2.*H  G.T.  HMAX,  STOP  DOUBLING  TEST 

9011  IF(NDUB)980. 9020.9020 

9020  IF (NVR) 980 .980 ,9021 

9021  IF(RAT-1. 19022, 9022, 980 

9022  IF(KS)954, 954, 9023 

9023  DO  903  1*1, IB 

IF (  ABS(DLT2(1)/AMAX1(  ABS (Y ( I ) ) , YHIN) ) »GE»VMIN)  00  TO  980 
C  IF(DA8S(DLT2(I)/DMAX1(0ABS(Y(I)),YMIN)).GE.VMIN)  GO  TO  980 

903  CONTINUE 

C  TEST  FURTHER  FOR  DOUBLING. 

DO  904  1=1 » IB 

IF  (  ABS(0LT1  (D/AMAX1  (  ABS  (Y  ( I) )  ,YMIN)  )  .GE.VMlN)  GO  TO  980 
C  IF(DABS(DLT1(I)/DMAX1(DABS(Y(I)) iYMIN) ) .GE.VHIN)  GO  TO  980 

904  CONTINUE 

C'  IF  HERE.  CURRENT  STEP  O.K..  FOR  DOUBLING. 

950  JMP4=H+H 

IF  (  A8SCTMP4)  .GT.HMAX)  GO  TO  3000 
C  IF (DABS (TMP4) .GT.HMAX)  GO  TO  3000  -  ~  ‘ 


DEQ40337 
DEQ40338 
DEQ40339 
DE040340 
DEQ40341 
”  ‘DEQ40342- 
DEQ40343 
DEQ40344 
DEQ40345 
OEQ40346 
DEQ40347 
DEQ40348 
DEQ40349 
DEQ40350 
DEQ40351 
DEQ40352 
DEQ40353 
DEQ40354 
DEQ40355 
DEQ40356 
DEQ40357 
DEQ40358. 
DEQ4.0359. 
DEQ40360 
DEQ40361 
DEQ40362 
DEQ40363 
DEQ40364 
DEQ40365 
DEQ40366 
DEQ40367 
DEQ4036Q 
DEQ40369 
DEQ40370 
DEQ40371 
DEQ40372 
DEQ40373 
DEQ40374 
DEQ40375 
DEQ40376 
DEQ40377 
DEQ40378 
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C  KJ*1  (SAVE  ALL  EQUATIONS  FOR  DOUBLING) 

KJ=  1 

952  KA=KC 
KC=KC+1 

GO  TO  (2010, 2020*2030, 2040,2050, 2060*2o70, 2080,2090) *KA 
954  KS=-KS 

GO  TO  1000 

980  IF(IB-Nl)981, 3002, 981 

981  IF(RXH-l.)3002,30n2,990 

990  IF (KS) 954,954,992 
992  KJsIB+1 

var=varx*vmin/rat 

DO  991  i«:<J,Nl 

IF (  ABS(DLT2(I)/AMAX1(  ABS ( Y ( I ) ) , yMIN) ) .GE. VAR)  GO  TO  3002 

991  CONTINUE 
NVR=-l 

GO  TO  952 

C  HALVE  H*  RETURN  Tn  R.K. 

1100  R=IR 

H=H/ (2«#R) 

IF (  ABS(H) .LT.HMIM)  GO  TO  3050 

C  IF (DARS(H) .LT.HHIm)  GO  TO  3050 

IDH=2 

rxh=rxa 

TST=*1. 

TEST=3. 

IF(NJH3) 1104, 1101,1101 

1101  IF (NSTRT) 11 02, 1108,1 106 

1102  TEST=2. 

C  RESTORE  INITIAL  CONDITIONS  AND  RESTART.  <TEST=2.) 

1104  T=TZ 
tst=test 

DO  1105  1=1, N1 
Y  ( I ) -T5 ( I ) 

YP ( I ) =T6 ( I ) 

1105  Y2P(I)=T7(I) 

1106  IKR=*1 
HA=H 

VMIN=ER/(H«H*100.) 
v‘max=vm  In*  1 00000  • 

GO  TO  502 

C  KA=1  SAVE  Y2P(0)  IN  T2  (ACCUM.  F1P(2H),  F2P(2H)) 


DEQ40379 

DEQ40380 

DEQ40381 

DEQ40382 

DEQ40383 

OEQ4O304 

PEQ40385 

DEG40386 

DEQ40387 

DEQ40388 

DEQ40389 

DEQ40390 

DEQ40391 

DE040392 

DEQ40393 

DEO40394 

DEQ40395 

DE040396 

DEO40397 

DEQ40398 

DEQ40399 

DEQ4O4OO 

DEO40401 

DEQ40402 

OE040403 

DE040404 

DEQ40405 

DEQ40406 

DEQ40407 

DEQ40408 

DEQ40409 

DEQ4O4I 0 

DE04041 1 

DEQ40412 

DEQ40413 

DEQ40414 

OEQ40415 

DEQ40416 

DE040417 

OE040418 

DEQ40419 

DEQ40420 
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2010  DO  2011  I»KJtNl  DEQ40421 

T2  < I ) *Y2P  < I )  DEQ40422 

T6U)sW0*Y2P(I)  DEQ40423 

2011  T7  (I1»V0*Y_2P(I>  OEQ40424 

<50  10  954  .  OEQ40425 

C5-5  -Kl*2-  SAVE" Y-2P i 2)-“IN  T3 *'  --fACCUM, -TlPTZH)-f-F2P^21flF)  --  -  --  DEQ40426 

2020  DO  2021  I*KJ»N1  DEQ40427 

T3 { I ) *Y2P ( I )  OEQ40428 

T 6;(.I.) »T6  ( I )  * W  1*Y2p ( I )  .  .r.-.  .  DEQ40429 

202 1  T71I )  «T7(  I )  f  V.l*Ygp  ( I )  .  DEQ40430 

GQl¥P^S4  :  ISg#;  '  sv  ...  DEQ40431 

C  Kl^‘3^ -5A\^-  Y2PW)^IN  T‘4  ■•TACCUMV‘TfPf2fl>*  F2Ft2HTr  -  ■ -v  DEQ40432 

2030  DO  2031  I=KJtNl  DEQ40433 

T4(I)sY2P<D  DEQ40434 

T6  M )  »T6  ( I  >'*W2*Y2»;(  I )  ---  --  DEQ40435- 

2031  "TT(I)»T7  (I)  +V2*Y2P.'(I)  <■'?->".  V.  DEQ40436 

GO.  TO  954-  .  ^  _  DEQ40437 

C  Kl*4  SAVE1':Y2P(fr)  IN  T5  <  ACCuMV'TIP^H)  *  F2PT2H) )  'ir-  -  DEQ40438 

2040  DO  2041  I«KJ»Nl  DEQ40439 

T5 ( I ) *Y2P ( I )  DEQ40440 

16  ( I  )  »T6  (.1 )  +W3#Y2P.(  I )  -  -  -  OEQ40441 

2041  17  (1 )  »T7  ( I  )4V3*Y2P  ( I )  ~  *  DEQ40442. 

GO  .  TO.  95.4* J  .  --  DEQ40443 

C-  Kl*5  ■  WCCUM.  F-1P(2H>  t  F2P(2H)  FOR  DOUBLING)  '  ‘  '  *  OEQ40444 

2050  H1*1,/(H*H)  DEQ40445 

H2*H1*H1  DEQ40446 

H3*1./RA-T  DEQ40447 

H4*H3*H3  DEQ40448 

DO  2051  ..  DEQ40449 

.  IF(I-IC>2053»2052V2053  DEQ40450 

2052  HlaHl*H3  DEQ40451 

H2*H2*H4  DEQ40452 

2053  T6 ( I ) *T6 ( I ) *W4*Y2P ( I ) +  YPT I ) *Hl  DEQ40453 

2051  T7(I)aT7(I>*V4*Y2p(I)+5**YP(l)*Hl*Y(I)*H2  DEQ40454 

GO  TO  954  DEQ40455 

C  Kl*6  (ACCUM.  F1P (2H) ♦  F2P(2H)  FOR  DOUBLING)  DEQ4045G 

2060  DO  2061  IsKJtNl  DEQ40457 

T6(I)sT6(I)+W5*Y2P(I)  DEQ40458 

2061  T7<I)sT7(I)+V5*Y2P(D  DEQ40459 

GO  TO  954  -  -  DEQ40460 

C  KJ*7  .(ACCUM.  F0JP(2H>,  E2P  C2H)-  :E0R§DpJjBblN6)-  -  .  DEQ40461 

2070  DO  2071  ^KJVNl  -  ‘  -5 -*S3* ^  --------  DEQ&0462 
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T6(I)=T6(I)+wo*Y?P{I) 

2071  T7(I)=T7(I)+V6*Y2P(I) 

GO  TO  954 

C  Kl=8  (ACCUM.  F1P (2H)  ♦  F2P(2H)  FOR 

2080  DO  2081  I=KJ*Nl 

T6<I)aT6(I)+W7*Y2P(I) 

208}  T7(I)=T7(I)+V7*Y2P(I) 

GO  TO  954 

C  Kl=9  (ACCUM.  F1P(2H),  F2P(2H)  FOR 
2090  DO  2091  I=KJ»N1 


IF(NH2)2095»3004»3004 
3004  IF ( JHH-3) 2095»300a»2095 
3006  NJH3=-1 
TZ=T 

DO  3010  1  =  1 t N 1 


DEQ40463 

DEQ40464 

DEQ4O465 

DOUBLING)  DEQ40466 

DEQ40467 

DEQ40468 

DEQ40469 

0EQ40470 

DOURLlNG)  DEQ40471 

DEQ40472 
DEQ40473 
0EQ40474 
DEQ40475 
DEQ40476 
DE040477 
DE040478 
DEQ40479 
0EQ40480 
DEQ40481 
DEQ40482 
DEQ40483 
DEQ40484 
DEQ40485 
DEQ4O486 
DEC40487 
OEQ40488 
OE040489 
DE040490 
DEQ40491 
DEQ40492 
DEQ40493 
DEQ40494 
DEQ40495 
DEQ40496 
DEQ40497 
DEQ40498 
DEQ40499 
OEQ40500 
PEQ40501 
DEQ40502 
DEO40503 
DEQ40504 


FlP(I)=T6(I>*w8*Y?P(I) 

F2P(I)sT7(I>*V8*Y?P<I) 

C  DOUBLE  H  IN  COWELL  MODE 

DLT8(I)=  2.*0LT6(I)-DLT7(I) 

DLT7(I)=  4.*DLT7<T)-<4.*DLT6(I)-DLT5(I) ) 

DLT6 ( I ) =  8.*DLT6 ( I ) - ( 12. *DL 15 ( I ) -6. *DLT4 ( I ) *DlT3 ( I ) ) 
DLT5(I)=16.*DuT5(T)-(32.*DLT4(I)-24.*DLT3(T)^8.*DLT2(T)-DLTl (I) ) 
DLT4(I)=  5.*ULT5(I)-(10,*DLT6(I)-108*DLT7(I)+B.*DlT8(T) ) 

1 ♦ ( Y2P ( I ) -T5 ( I ) ) 

DLT3(I)=  6.*DLT4(I)-<l5.*OLT5<I)-20.*DLT6(I)M5.*nLT7(l) 
1-6.*DLT«(I) )-(Y2P(I)-T4(I) ) 

DLT2(I)s  7.*ULT3U)-<21.*0LI4<I)-35.*DLT5(I)+35.*DLT6(I> 
1-21#*0LT7<D +/.*DLT8(I) ) + (Y2P( I ) -T3 ( I) ) 

DLT1 (I)=  8.#DLT2(I)-(28.*DLT3(I)-56.*DLT4(I) ♦70.*DLT5(I) 
1-56.*DLT6 ( I ) +28.*DLT7 ( I ) -8.*0LT8( I) ) - ( Y2P ( T ) -T2 ( I ) ) 

2091  CONTINUE 

C  END  OF  DOUBLING 

IF(NVR)2093t2u92t2092 

2092  IDH=3 
VMIN=VMlN/4. 

VMAX=VMAX/4. 

HsH+H 

HA=H 
NH2=NHH 
GO  TO  2095 
3000  NDUBs-1 
3002  NJH3S+1 
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T5 ( I ) =Y { I ) 

DEQ40505 

T6(I)=YP(I) 

DEQ40506 

3010 

T7  < I ) sY2p  ( I ) 

DEQ40507 

GO  TO  2095 

DEQ40508 

3050 

NH2=-1 

DEQ40509 

IF (NSTRT) 3052* 1000 *1000 

DEQ40510 

3052 

TEST=13. 

DEQ40511 

GO  TO  261 

DEQ40512 

2093 

RAT=2.«RAT 

DEQ40513 

K=K*K 

OEQ40514 

RKH=kKH-l. 

DEQ40515 

I0H=4 

DEQ40516 

2095 

KC=1 

DEQ40517 

KS=  1 

DEQ40518 

C 

STEP-SIZE  O.K.t  PROCEED  WITH  COWELL  STEP. 

DEQ40519 

1000 

TEST=-1. 

DEQ40520 

H=HA 

DEQ40521 

NSTRT=.l 

DEQA0522 

HH=H*H 

DEQ40523 

NX=  IB 

DEQ40524 

C 

complete  k  steps  for  first  ib  eons. 

DEQ40525 

DO  1022  K1=1*K 

DEQ40526 

T=T  *H 

DEQ40527 

IF  (K-Kl ) 6660 .6660 ,6661 

DEQ40526 

6660 

MX=N1 

DEQ40529 

TEST=*1. 

QEQ40530 

6661 

DO  1001  1  =  1 » NX 

DE040531 

IF  (I-IC) 6666.6665,6666 

DEQ40532 

6665 

h=ha*rat 

DEQ40533 

HH=H*H 

0EQ40534 

6666 

T1  ( I ) =Y2P  < I ) 

DEQ40535 

1001 

Y ( I >  =HH* (F2P ( I ) ♦£0*Y2P ( I ) ♦Al*OLT8 < T  > ♦ A2*DLT7 ( I >  + A3*DLT6 ( f ) 

DEQ40536 

1  ♦A4*DLT5 ( I > ♦A5*0L  T4 (I ) *A6*OLT3 (I ) +A7*DLT? (I) ♦A8*DLT1 ( I )  ) 

DEQ40537 

H=HA 

DEQ40538 

IF  C 1 1 )  1007.1005.1005 

DEQ40539 

1005 

DO  1006  1=1 .NX 

DEQ40540 

IF(I-IC) 1006.6667,1006 

DE040541 

6667 

h=ha*rat 

0EO40542 

1006 

YP<I>=H*<FlPU)  ♦Ap0*Y2P<I)*AP1*DLT8(I)*AP2*DLT7<I>+AP3«DLT6(I> 

DEQ40543 

1  ♦AP4*DLT5n>*AP5*DLT4(I>*AP6*DLT3<I>*AP7*DLT2<lUAP8*f)LTUin 

DEQ40544 

1007 

ASSIGN  1015  TO  IFLAG 

DEQ40545 

IA«-1 

DEQ40546 
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1010  CALL  DAUX 
H=HA 
HH=H»H 

1009  DO  1013  I s:  1 1 N X 

IF  { J  —  I  C )  6B&9  » 0568  « 6669 
6668  h=HA»RAI 
HH=H#H 

b669  TMP1  =  V2P(I)”IJ  '•  I ) 

TMP2=TMPl-0LiH(l) 

7MP3-TMHH-DL r / ( I ) 

TmPX=TmP3-DLIo( () 

TMPB-TMP4-0I.  I  J  ( I ) 

TMP6=TMP5-0Lr4(l) 

TMP7=TMP6-0Lfj(l) 

TMP8=TMP7-l)Ll  2  (  ( ) 

Y(I)=HH»(F3P( 1) ♦  «f<#Y2P(I) *bl*TMPl*8?*rMP2*R3*TMP3*B4*TMP4*R5#TMP5 
1  +B6>#TMP6  +  H7ttlMP/  +  bH*T'-1PB) 

1011  YP!I>  =H#  (FlP  1  a  )  ♦BP0*Y2P(I)*bPl*TMPUHP?*TMP2  +  RP3ttrMP3*RP4*TMP4 
1  +BPS*TMPR+HPo«  rMP6  +  BP7*TMP7  +  RP8»TMPd> 

1013  CONTI N" It 

00  TO  IFLAOo (1018.1020) 

1015  ASSIGN  1020  Tu  IFi  AG 
IA=+l 

GO  TO  1010 

1020  DO  1021  I s 1 » NX 
TMPl=Y2P(I)~ri  (I) 

TMP2=TMPl-l)L  I  8  ( I ) 

DLTB ( I ) STMP  J 
TMPl  =  TMP2-m.r  Ml) 

DLT7(I)=IMP2 

TMP2=TMPi-OLTu(I) 

0LT6( I ) =  T  MP 1 
TMPl=TMP2-0Lrb(I) 

DLT5 ( r)=lMP2 
ImP?=TNP1-DL  T'h  (I ) 

DLT4 ( T ) = l MP1 
TMP1=TMP2-0LTJ ( t ) 

DLT3<I)-"TmP2 

DLT 1 (I)=TMP)-OLr?(I) 

OLT2( I >  =  TMP 1 
F1P(I)=F1P( U  *722(1) 

1021  F2P { I i =F2P ( I ) *F IP  < 1 ) 


DEQ40547 

DEQa.0548 

DEQ40549 

DE040550 

OEQ40551 

UEQ40552 

DEQ40553 

DE040554 

UEQ4055R 

OEQ40556 

DEO40557 

DEO40558 

DE040559 

DEO405C0 

DEQ40561 

DE040562 

OEQ4O563 

DEQ40564 

DEQ4O565 
DEO40566 
OEQ40567 
DEQ4O560 
DE04C569 
DEQ4057O 
0EQ4057 1 
OEQ40572 
DEQ40573 
DEQ40574 
OEQ4057S 
DEO40576 
OEQ40577 
DEQ40578 
DEQ40579 
DEQ40580 
DEQ40581 
DE040582 
DEQ40583 
DEO40>384 
DEQ40585 
DEQ40586 
DEQ4O587 
OEQ40588 
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DEQ40589 

TEST=TFST-1 •  UEQ40590 

K=HA  DEQ40591 

HH=HA#HA  OEQ40592 

1022  CONTINUE  DE040593 

TEST=-1 •  OEQ40594 

GO  TO  201  DEQ40595 

END 
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